Characterizing mitochondrial features in osteoarthritis through integrative multi-omics and machine learning analysis

Purpose Osteoarthritis (OA) stands as the most prevalent joint disorder. Mitochondrial dysfunction has been linked to the pathogenesis of OA. The main goal of this study is to uncover the pivotal role of mitochondria in the mechanisms driving OA development. Materials and methods We acquired seven bulk RNA-seq datasets from the Gene Expression Omnibus (GEO) database and examined the expression levels of differentially expressed genes related to mitochondria in OA. We utilized single-sample gene set enrichment analysis (ssGSEA), gene set enrichment analysis (GSEA), and weighted gene co-expression network analysis (WGCNA) analyses to explore the functional mechanisms associated with these genes. Seven machine learning algorithms were utilized to identify hub mitochondria-related genes and develop a predictive model. Further analyses included pathway enrichment, immune infiltration, gene-disease relationships, and mRNA-miRNA network construction based on these hub mitochondria-related genes. genome-wide association studies (GWAS) analysis was performed using the Gene Atlas database. GSEA, gene set variation analysis (GSVA), protein pathway analysis, and WGCNA were employed to investigate relevant pathways in subtypes. The Harmonizome database was employed to analyze the expression of hub mitochondria-related genes across various human tissues. Single-cell data analysis was conducted to examine patterns of gene expression distribution and pseudo-temporal changes. Additionally, The real-time polymerase chain reaction (RT-PCR) was used to validate the expression of these hub mitochondria-related genes. Results In OA, the mitochondria-related pathway was significantly activated. Nine hub mitochondria-related genes (SIRT4, DNAJC15, NFS1, FKBP8, SLC25A37, CARS2, MTHFD2, ETFDH, and PDK4) were identified. They constructed predictive models with good ability to predict OA. These genes are primarily associated with macrophages. Unsupervised consensus clustering identified two mitochondria-associated isoforms that are primarily associated with metabolism. Single-cell analysis showed that they were all expressed in single cells and varied with cell differentiation. RT-PCR showed that they were all significantly expressed in OA. Conclusion SIRT4, DNAJC15, NFS1, FKBP8, SLC25A37, CARS2, MTHFD2, ETFDH, and PDK4 are potential mitochondrial target genes for studying OA. The classification of mitochondria-associated isoforms could help to personalize treatment for OA patients.


Introduction
Osteoarthritis (OA) is a chronic and disabling joint disease characterized predominantly by the progressive destruction of articular cartilage.It is the leading cause of declining quality of life among elderly individuals worldwide.According to the latest statistics from the World Health Organization, OA affects over 500 million people globally (1).The absence of a cure for OA can be attributed, in part, to an incomplete comprehension of the underlying pathological mechanisms governing its onset and progression.Consequently, it becomes imperative to enhance our understanding of the associated signaling pathways and critical molecules involved in OA development.This knowledge will facilitate the identification of therapeutic targets and aid in the design of effective treatments.
Mitochondria, as organelles within eukaryotic cells, produce adenosine triphosphate (ATP) through oxidative phosphorylation, thereby serving as the cell's energy source (2).Additionally, they play a crucial role in cellular metabolism and maintaining cellular homeostasis.Impaired mitochondrial function reduces the capacity for oxidative phosphorylation, leading to an upsurge in reactive oxygen species (ROS) production.This accumulation of ROS can cause cellular damage or apoptosis (3).Mitochondria, as essential organelles in chondrocytes, play critical roles in cellular metabolism, proliferation, and apoptosis.Growing evidence suggests that mitochondrial dysfunction and disrupted energy metabolism are closely associated with the onset and progression of OA.However, the exact mechanisms remain unclear at present.
Research indicates that impaired mitochondrial function leads to increased levels of ROS, resulting in chondrocyte apoptosis and promoting the development of OA (4).It is noteworthy that oxidative stress is a key factor in causing mitochondrial DNA damage, impairing mitochondrial respiratory function, and activating mitochondria-mediated cell death pathways (5).Inflammatory cytokines like IL-1b and tumor necrosis factoralpha (TNF-a) have been reported to reduce mitochondrial activity and ATP production, impairing mitochondrial respiration in chondrocytes and contributing to mitochondrial dysfunction in OA (6).Additionally, the balance between mitochondrial fission and fusion is essential for preserving proper mitochondrial function.Increased fission and decreased fusion can cause mitochondrial fragmentation and functional impairment, leading to reduced ATP production and increased ROS generation (7).Disrupted mitochondrial segregation during mitosis has also been implicated in mitochondrial dysfunction in OA (4).Additionally, alterations in mitochondrial metabolism may disrupt the cellular redox balance and contribute to the accumulation of ROS in OA.During the progression of OA, chondrocytes and synovial cells tend to shift their mitochondrial metabolism from oxidative phosphorylation to glycolysis, primarily regulated by the AMPK-MTOR pathway (8).Increasing evidence suggests that the loss of mitochondrial quality control homeostasis contributes to cartilage damage in the onset and progression of OA (9).Mitochondrial quality control is a crucial mechanism that coordinates various mitochondrial biological functions.Autophagy, a vital cellular mechanism for maintaining homeostasis, plays a critical role in chondrocytes and is closely linked to the pathogenesis of OA.Mitophagy, specifically, plays a crucial role in degrading, clearing, and recycling dysfunctional mitochondria, thus controlling mitochondrial quality and maintaining cellular homeostasis (10).Animal studies have confirmed that inhibiting mitophagy promotes cartilage degradation (11).Emerging evidence suggests that mitophagy plays a significant regulatory role in skeletal disorders, suggesting that modulating mitophagy levels could be a novel strategy for treating skeletal-related diseases (12).In summary, multiple factors, including oxidative stress, exacerbated inflammation, mitochondrial fission and fusion, disrupted mitosis, metabolic dysregulation in chondrocytes, and impaired autophagy, contribute to mitochondrial dysfunction in OA.Therefore, comprehensive understanding of the mechanisms by which mitochondrial dysfunction leads to the onset and progression of OA is crucial for developing effective strategies for its treatment.
By integrating bulk RNA, weighted gene co-expression network analysis (WGCNA), immune infiltration analysis, and molecular subtyping, we have investigated the potential alterations in mitochondria-related genes, their impact on pathways, and the role of immune cells.Additionally, we have identified targeted therapies, immune microenvironment, and pathway heterogeneity within different mitochondrial metabolic subgroups.Furthermore, we have developed and validated a predictive model for OA using machine learning and deep learning.Single-cell analysis revealed single-cell expression distribution patterns and putative temporal changes of hub mitochondria-related genes.Finally, real-time polymerase chain reaction (RT-PCR) verified the expression of hub mitochondria-related genes.

Data acquisition and pre-processing
We retrieved OA-related chip datasets from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), which included GSE117999, GSE51588, GSE55235, GSE55457, GSE57218, GSE82107, and GSE98918.Mitochondria-related genes were obtained from the MitoCarta3.0 database.Overview of the information on analyzed datasets is provided in Supplementary Table 1.During the processing of the data, we excluded ineligible samples.During the processing of the chip datasets, we matched probe and gene names based on the annotation information of each GPL platform.In cases where multiple probes corresponded to the same gene, we selected the probe with the highest expression level for retention.Following this, we normalized the expression matrix using the `normalizeBetweenArrays`function and performed log2 transformation for datasets that required it.Next, we extracted the common genes shared by all seven datasets.To address batch effects and platform variations, we employed the ComBat method from the R package "sva" (13,14) to normalize expression values across different batches or platforms.The success of batch effect removal was evaluated through principal component analysis (PCA).

Differential expression analysis and enrichment analysis
Differential gene expression analysis was performed using the limma R package (15,16) X to compare the expression profiles between OA samples and normal samples.The criteria for selecting differentially expressed genes (DEGs) associated with mitochondria were set as a p-value < 0.05.Next, we evaluated the activity scores of each pathway using the single-sample gene set enrichment analysis (ssGSEA) (17) and gene set enrichment analysis (GSEA) (18).Then, we utilized the clusterprofiler R package (19,20) to perform functional and pathway analysis, including Gene Ontology (GO) (21) and Kyoto Encyclopedia of Genes and Genomes (KEGG), for the DEGs associated with mitochondria.

WGCNA analysis
In the WGCNA analysis, we used the mitochondria-related DEGs as the background gene set.We calculated the Gene Set Variation Analysis (GSVA) (22) scores for each sample and used these GSVA scores as traits for the WGCNA analysis.On the basis of selecting all genes, the determination of the maximum R2 (power = 4) was conducted using the soft threshold of an unweighted network.Each module is required to contain at least 30 genes, and the topological overlap matrix similarity was utilized to evaluate the distance between gene pairs.Hierarchical clustering analysis was performed using the average linkage method and dynamic tree cut algorithm to construct a clustering tree and divide the genes into different modules.Next, we utilized the clusterprofiler R package (23) to perform GO and KEGG analysis, on the modules that were most correlated with the GSVA scores.String database to construct protein interaction networks for modular genes.

Identification of hub mitochondriaassociated DEGs
To identify hub genes among the mitochondria-related DEGs, we employed the Lasso regression, ridge regression, and elastic net regression for feature selection.Additionally, we utilized various methods such as SVM-RFE, Random Forest, Gradient Boosting Machine (GBM), and XGBoost to rank the importance of genes.We selected the top 30 genes based on their importance scores derived from these methods.The criteria for identifying hub genes among the mitochondria-related DEGs were set as the genes that appeared consistently across all seven algorithms mentioned above.

Predictive model and alignment diagram
To determine the best machine learning model for predicting OA, we selected the expression profiles of hub mitochondriaassociated DEGs as input variables, and built seven machine learning models including logistic regression, SVM, kknn, random forest, LDA, Naive Bayes, and decision tree for prediction using the "mlr3verse" R package.The reliability of the models was also verified using bootstrap method.The genes were then evaluated using multifactorial logistic regression analysis, and the area under the receiver operator characteristic (ROC) curve (AUC) was calculated using the "ROCR" package to assess the value of these genes in OA diagnosis.We used a nomogram to predict the likelihood of OA, and plotted correction curves and decision curves to analyze the stability of the model.

Deep learning
Deep learning is a type of neural network that has more hidden layers and can understand more complex information structures in the data.Multi-Layer Perceptron (MLP) is a common form of neural network, also known as feedforward neural network.We compute the confusion matrix and visualize the MLP diagnostic results, extract the weights and biases of each layer, and finally generate the network structure and visualization.

GSEA and ssGSEA analysis of hub mitochondria-associated DEGs
Using the mitochondrial gene set and KEGG pathway gene set as the background gene set, we performed GSEA analysis of hub mitochondria-related DEGs using the clusterprofiler R package.The ssGSEA algorithm was also utilized to calculate the score of the hallmark pathway set, and the correlation between the hub mitochondria-related DEGs and pathways was calculated using the spearman algorithm.

Immune infiltration analysis
The cibersort (24) algorithm was utilized to calculate the proportion of immune cells in OA, and subsequently, the spearman algorithm was employed to assess the correlation between hub mitochondria-associated DEGs and immune cells.Additionally, the spearman algorithm was used to determine the correlation between hub mitochondria-associated DEGs and immune-related genes.

Constructing gene-disease networks and mRNA-miRNA networks
The DisGeNET R package (25) was utilized to analyze diseases associated between genes, while multiple miRAN databases were integrated to predict the miRAN of genes.

Genome−wide association study
The GeneAtlas database (http://geneatlas.roslin.ed.ac.uk/) serves as a comprehensive resource that utilizes data from the UK Biobank cohort to provide extensive information on associations between hundreds of traits and millions of genetic variations.The database contains information from the UK Biobank, including data from 452,264 individuals.It covers 778 phenotypes and encompasses a vast number of genetic loci, totaling approximately 30 million.It offers valuable genetic and phenotypic data to researchers, facilitating a deeper understanding of the relationships between genes and traits and advancing research in related fields.

Consensus clustering
Consensus clustering (26, 27) is a comprehensive clustering method designed to identify subgroups within a dataset.In this study, we performed clustering analysis using the K-means algorithm on the expression profiles of mitochondria-related DEGs in OA samples.We specified a maximum of 10 clusters and determined the final number of clusters based on the consistency matrix and cluster consensus score, which exceeded 0.8.Subsequently, we further analyzed the differential expression levels of the mitochondria-related DEGs and immune cells across different subgroups.

Drug analysis
Utilizing the results from the subtype differential expression as disease indicators and the drug data in the CMap database as drug markers, we implemented the eXtreme Sum (XSum) methodology for feature alignment.This was done to identify drugs that are suitable for particular subtypes.

GSEA, GSVA and protein function analysis of subtypes
To look for subtype-related fe and dynamic tree-cutting methods.This approach facilitated the creation of a dendrogram and the assignment of genes to their respective modules.After considering their similarities, we merged the original 11 modules and ultimately identified six functional modules.To further analyze the key modules, we utilized the ClusterProfiler R package to perform GO and KEGG functional analyses.

Analysis of expression levels in normal human tissues
Expression of hub mitochondria-related genes in human tissues was analyzed using the harmonizome database (https:// maayanlab.cloud/Harmonizome/).

Single-cell analysis
We used the GSE133449 dataset for single-cell level expression distribution and pseudotime analysis.We processed the single-cell RNA sequencing data using the Seurat R package (28).Cells expressing more than 200 but not exceeding 2,500 genes were selected.The "FindVariableGenes" function was employed to identify highly variable genes, followed by PCA.For dimensionality reduction and visualization of the single-cell data, we utilized Uniform Manifold Approximation and Projection (UMAP).To generate a plot of the single-cell data, we used the "DimPlot" function, and for visualizing gene expression patterns, we employed the "FeaturePlot" function.Cell type annotation was performed using the 'HumanPrimaryCellAtlasData()' function from the 'celldex' package.Subsequently, we calculated senescence scores for each cell type and evaluated differences in senescence scores between the normal and OA groups for the same cell types.Further visualization of these differences can be achieved using the "FeaturePlot" function.Note: 'HumanPrimaryCellAtlasData()' is not actually a function in the 'celldex' package.It is used here as an example for cell type annotation.Please refer to the appropriate package or resource for cell type annotation in your specific analysis.

Pseudotime analysis
We extracted a subset of chondrocyte cells from the single-cell data for pseudotime analysis.Initially, we performed redimensionality reduction and clustering specifically on the chondrocyte cells.Next, we utilized the Monocle R package for pseudotime analysis.We selected cells for subsequent pseudotime analysis based on average expression greater than 0.1 and empirical dispersion greater than one-fold of the fitted dispersion.To reduce dimensionality and order the cells in the pseudotime analysis, we used the 'DDRTree' method from the 'reduceDimension' function.Additionally, We examined the expression variations of mitochondria-related genes across distinct clusters during the cell differentiation trajectory.

RT-PCR validation
Human chondrocyte cells were obtained from Wuhan Saos Technology Co., Ltd.These cells were cultured in DMEM/F12 medium supplemented with 10% fetal bovine serum.To induce inflammation, the model group's cells were treated with a culture medium containing IL-1b at a concentration of 10 ng/ml for 24 hours.Total RNA extraction was performed using the QIAzol reagent kit, followed by cDNA synthesis through reverse transcription using oligo-dT primers.Subsequently, gene amplification of the cDNA was carried out using an RT-PCR machine, following an initial denaturation at 95°C for 5 minutes, followed by 40 cycles of denaturation at 95°C for 1 minute, annealing at the optimal temperature of 60°C for 30 seconds, and extension at 72°C for 1 minute.The relative expression level of each specific gene was calculated using the 2-DD method.Please refer to Table 1 for the primer sequences of the human genes utilized.

Data processing and enrichment analysis
We integrated and processed data from GSE117999, GSE51588, GSE55235, GSE55457, GSE57218, GSE82107, and GSE98918, resulting in 66 control samples and 125 OA samples, encompassing information on 10827 genes.As depicted in Figure 1A, batch effects were observed across the seven different datasets.However, after eliminating these batch effects, consistent results were obtained (Figure 1B).This indicates that we successfully eliminated batch effects through cross-platform standardization methods.Differential expression analysis revealed 245 mitochondria-related DEGs, including 72 up-regulated and 173 down-regulated mitochondria-related DEGs (Figure 1C).The ssGSEA results demonstrated a significant decrease in mitochondrial activity in OA (Figure 1D), which was further validated by GSEA analysis (Figure 1E).Enrichment analysis revealed significant enrichment of mitochondria-related DEGs in processes such as ribose phosphate metabolic process, mitochondrial gene expression, mitochondrial transport, Pathways of neurodegeneration -multiple diseases, Parkinson disease, and Prion disease (Figure 1F).WGCNA analysis indicated a significant correlation between the black module and mitochondrial scores (Figure 1G). Figure 1H displays the Module Membership (MM) and Gene Significance (GS) for each gene in the black module.Both BP and KEGG analyses showed significant enrichment of genes in the black module in processes such as cytokine-mediated signaling pathway, positive regulation of cell adhesion, response to xenobiotic stimulus, Cytokine-cytokine receptor interaction, Calcium signaling pathway, and Ras signaling pathway (Figure 1I).The protein-protein interaction network of genes in the black module is illustrated in Figure 1J.

Machine learning identification of hub mitochondria-related DEGs
We successfully utilized lasso regression (Figure 2A), ridge regression (Figure 2B), and elastic regression (Figure 2C) to screen for important mitochondria-related DEGs.Additionally, we employed Support Vector Machine -Recursive Feature Elimination (SVM-RFE) to select the top 30 important mitochondriarelated DEGs (Figure 2D).Moreover, we used random forest

Primer
Sequence  2E, F), Xgboost (Figure 2G), and GBM (Figure 2H) to independently identify the top 30 mitochondria-related DEGs.Subsequently, we obtained the intersection of these seven machine learning results, resulting in nine hub mitochondria-related DEGs, namely SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 (Figure 2I).The heatmap displayed the expression patterns of SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 in OA and control samples (Figure 2J).Box plots illustrated the expression levels of these genes in OA and control samples (Figure 2K).The correlation and protein-protein interaction information of these genes are presented in Figures 2L and M, respectively.Functional similarity analysis revealed that PDK1 holds a significant position among these genes (Figure 2N).The chromosomal locations of these genes are shown in Figure 2O.

Column chart and neural network model construction
After comparing seven machine learning models, we found that the random forest model exhibited the highest AUC value (Figure 3A).Additionally, the random forest model demonstrated good sensitivity and specificity (Figure 3B).The AUC of the random forest model was 0.968 (Figure 3C), indicating its ability to accurately identify OA patients.Internal validation further confirmed the reliability of the model (Figure 3D).In the column chart, each feature variable corresponds to a specific score, and the sum of all feature scores represents the probability of having OA (Figure 3E).The calibration curve validated the accuracy of the column chart in diagnosing OA (Figure 3F).Decision Curve Analysis (DCA) showed that the column chart provided certain benefits for OA patients in clinical applications (Figure 3G).During the construction of the neural network model, Figure 3H displayed the accuracy and loss rate for each of the 100 training iterations.As indicated by the Kappa test (Figure 3I), the confusion matrix exhibited good consistency.Figure 3J presented a comparison between the predicted and actual results.The structure of the neural network model is depicted in Figure 3K.The model achieved an AUC of 0.961, demonstrating high diagnostic performance (Figure 3L).

Enrichment analysis of hub mitochondriarelated DEGs
Except for SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 were significantly activated in mitochondrial pathways (Figure 4A). Figure 4B provides detailed information on the top six pathways that showed significant enrichment when using the KEGG pathway gene set as the background gene set for hub mitochondria-related DEGs. Figure 4C illustrates the correlation between hub mitochondriarelated DEGs and hallmark pathway gene sets.

Immune infiltration analysis
The heatmap in Figure 5A depicts the correlation between hub mitochondria-related DEGs and immune cells.Subsequent selection of results with significant correlations is presented in the scatter plot, providing a detailed visualization of the correlation between each hub mitochondria-related DEG and individual immune cells (Figure 5B).SIRT4 exhibits positive correlation with Dendritic cells resting and T cells CD4 naive.DNAJC15 displays positive correlation with Macrophages M2, NK cells activated, and T cells gamma delta, while exhibiting negative correlation with NK cells resting and T cells regulatory (Tregs).NFS1 demonstrates positive correlation with Plasma cells and negative correlation with Mast cells resting and Dendritic cells resting.FKBP8 exhibits negative correlation with Dendritic cells resting and T cells gamma delta.SLC25A37 shows positive correlation with Neutrophils, Monocytes, NK cells resting, and Plasma cells, while demonstrating negative correlation with Dendritic cells resting, Macrophages M2, and T cells gamma delta.CARS2 exhibits negative correlation with T cells CD8.MTHFD2 shows negative correlation with NK cells resting.ETFDH displays positive correlation with Monocytes and NK cells resting, while showing negative correlation with Mast cells resting.PDK4 demonstrates negative correlation with Macrophages M1, Monocytes, and T cells CD8, while displaying negative correlation with Macrophages MO.Further analysis reveals the correlation between hub mitochondria-related DEGs and immunerelated genes, as depicted in Figure 5C.

Gene-disease network and mRNA-miRNA network
By analyzing the associations between genes and diseases, we constructed a gene-disease network (Figure 6A).For hub mitochondria-related DEGs, we selected five databases that encompassed the maximum number of miRNAs (Figure 6B) and successfully constructed an mRNA-miRNA network (Figure 6C).

Consensus clustering analysis
We performed classification of OA samples based on hub mitochondria-related DEGs.Using the consensus matrix as a similarity matrix, we determined the final subtypes.Based on the consensus clustering results, cumulative distribution function (CDF) plot, relative changes in CDF curves, and consensus clustering scores, we selected k=2 as the optimal value, dividing the OA samples into two distinct subtypes (Figures 8A-D).In subtype 1, CD8 T cells showed significantly lower expression, while M0 macrophages exhibited significantly higher expression (Figure 8E).Drug analysis revealed that STOCK1N.35696,fasudil, MK.886, and X4.5.dianilinophthalimide were among the top five potential drugs for treating subtype 1 patients (Figure 8F), whereas clofibrate, MS.275, NU.1025, imatinib, and butein were among the top five applicable drugs for treating subtype 2 patients (Figure 8G).Figures 8H and 8I respectively illustrate some biological processes and KEGG pathways involved in the subtypes.

Protein-level enrichment analysis and GSVA analysis of Subtypes
Protein pathway analysis revealed that subtype 1 was primarily enriched in the PI3K-Akt signaling pathway (Figure 9A), while subtype 2 showed enrichment in transcription factors (Figure 9B).GSVA analysis indicated that both subtype 1 and subtype 2 were associated with numerous metabolic functions or pathways in terms of biological processes (BP) (Figure 9C) and KEGG pathways (Figure 9D).

WGCNA analysis of subtypes
Using the information from subtype 1 and subtype 2, we conducted WGCNA analysis.The soft-thresholding power was set to 4 (Figure 10A), with a correlation coefficient threshold of 9.This resulted in the generation of 27 modules (Figure 10B).Among these modules, the blue module exhibited the strongest correlation with the subtypes (Figure 10C).Inter-module correlations are illustrated in Figure 10D.Module-trait analysis focused on the blue module, displaying the distribution of MM and GS for genes within this module (Figure 10E).Regarding BP, enrichment within the module primarily involved signal release, response to peptide hormone, and Enrichment analysis of hub mitochondria-associated DEGs.(A) GSEA analysis for SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 using the mitochondria gene set as the background gene set.(B) GSEA analysis for SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 using the KEGG pathway gene set as the background gene set.(C) ssGSEA analysis for SLC25A37, MTHFD2, SIRT4, DNAJC15, ETFDH, PDK4, CARS2, FKBP8, and NFS1 using the hallmark pathway gene set as the background gene set.
second-messenger-mediated signaling (Figure 10F).In terms of KEGG pathways, the module genes were predominantly enriched in Neuroactive ligand-receptor interaction, cAMP signaling pathway, and Cytokine-cytokine receptor interaction (Figure 10F).

Human tissue expression atlas
Utilizing the Harmonizome database, we analyzed the expression of hub mitochondria-related DEGs in human  11A-I),with a particular focus on bone and cartilage cells.

RT-PCR validation results
The OA group exhibited higher expression levels of DNAJC15, FKBP8, NFS1, and SIRT4 genes compared to the control group.Conversely, the OA group showed lower expression levels of CARS2, ETFDH, MTHFD2, PDK4, and SLC25A37 genes in comparison to the control group (Figure 13).

Discussion
OA is typically caused by abnormal metabolism and death of chondrocytes.The functionality of mitochondria is crucial for maintaining a stable intracellular environment in chondrocytes.Impaired mitochondrial function, a significant factor in cellular metabolic abnormalities, contributes to the onset and progression of OA.Recent studies have also reported that dysfunctional mitochondria in chondrocytes can lead to the generation of ROS, oxidative stress, inflammation, and cell apoptosis (29).Therefore, alterations in mitochondrial function are associated with the pathological processes of OA.
In this study, we systematically analyzed genes associated with mitochondrial metabolism in OA.Through differential expression analysis, we identified a set of differentially expressed genes linked to mitochondrial metabolism.Furthermore, employing the SsGSEA algorithm, we evaluated the expression levels of the mitochondrial metabolism-related gene set and observed significant downregulation in OA.This finding highlights the influence of mitochondrial metabolic dysfunction on the development of OA.Through enrichment analysis, we discovered that the main biological processes involving mitochondrial metabolism-related genes included mitochondrial gene expression, mitochondrial transport, mitochondrial translation, mitochondrial transmembrane transport, and mitochondrial membrane organization.Furthermore, we identified significant signaling pathways associated with mitochondrial metabolism-related genes, including cytokinemediated signaling pathway, positive regulation of cell adhesion, response to xenobiotic stimulus, cytokine-cytokine receptor interaction, calcium signaling pathway, Ras signaling pathway, cGMP-PKG signaling pathway, hematopoietic cell lineage, and HIF-1 signaling pathway.These findings serve as a theoretical foundation for further exploring the underlying mechanisms of mitochondrial metabolism in OA.
Through the use of seven machine learning algorithms, we identified nine mitochondrial metabolism-related genes, namely SIRT4, DNAJC15, NFS1, FKBP8, SLC25A37, CARS2, MTHFD2, ETFDH, and PDK4, which are closely associated with the occurrence of OA.Based on these nine genes, we constructed a prediction model for OA, which can be utilized for diagnosis and treatment purposes.Previous studies have confirmed that sirtuins (SIRT) play a crucial role in the regulation of cellular metabolism and are associated with various diseases.SIRT4, a member of the SIRT protein family, is specifically localized within mitochondria and maintains mitochondrial function and homeostasis (30).Overexpression of SIRT4 contributes to the inhibition of inflammation and oxidative stress in OA (31).SIRT4 shows a positive correlation with the transcriptional activity of hypoxiainducible factor-2a (HIF-2a) in chondrocytes.HIF-2a transcription upregulates NAMPT in articular chondrocytes, leading to the development of OA in mice.The interplay between HIF-2a and the NAMPT-NAD(+)-SIRT axis in chondrocytes is associated with HIF-2a or NAMPT-induced cartilage destruction in OA (32).Overexpression of SIRT4 in chondrocytes enhances mitochondrial autophagy, preventing cellular senescence and degradation (33).DNAJC15, also known as methylation-controlled J protein (MCJ), is an inner mitochondrial membrane protein serving as an endogenous inhibitor of respiratory chain complex I. Its depletion results in enhanced mitochondrial respiration and elevated reactive oxygen species levels (34).It is part of the mitochondrial protein transport mechanism and is involved in cell death processes (35).Lack of MCJ affects the pathological physiology resulting from alterations in mitochondrial metabolism (36).NFS1 (cysteine desulfurase) provides persulfide for the biosynthesis of ironsulfur (FeS) clusters in mitochondria (37).FeS clusters are essential for the activity of many cellular proteins.FKBP8 belongs to the FK506-binding protein (FKBP) family and is a membrane-associated protein typically found in mitochondria.FKBP8 regulates autophagy by interacting with the VPS34 lipid kinase complex (38).Solute carrier family 25 member 37 (SLC25A37), also known mitochondrial iron transporter-1, is a member of the solute carrier family located in the inner mitochondrial membrane (39).It mediates the uptake of mitochondrial iron.Persulfides are sulfur-containing organic compounds that possess diverse physiological functions.The synthesis of persulfides is primarily mediated by mitochondrial cysteinyl-tRNA synthetase (CARS2).Persulfides mediate mitochondrial biogenesis and bioenergetics.The function of mitochondria in biogenesis and bioenergetics is also supported and upregulated by persulfides derived from mitochondrial cysteinyl-tRNA synthetase (CARS, also known as CARS2) (40).Methylenetetrahydrofolate dehydrogenase 2 (MTHFD2) is a mitochondrial enzyme that participates in one-carbon metabolism.It plays a crucial role in maintaining the integrity of the mitochondrial respiratory chain and preventing mitochondrial dysfunction (41).The electron transfer flavoprotein (ETF), consisting of a and b subunits (ETF a and ETF b, respectively), functions as a mitochondrial matrix protein.It facilitates electron transfer between various mitochondrial dehydrogenases and the primary respiratory chain through the action of electron-transferring flavoprotein dehydrogenase (ETFDH) (42).Pyruvate dehydrogenase kinase 4 (PDK4) inhibits the development of OA through activation of the PPAR pathway (43).Overexpression of PDK4 is also sufficient to promote mitochondrial fission.The PDK4-SEPT2-DRP1 axis has been identified as a regulator of mitochondrial function, bridging the gap between cellular bioenergetics and mitochondrial dynamics (44).However, apart from SIRT4 and PDK4, no other studies have been found examining the remaining seven genes in relation to OA.Our study further confirms the association between mitochondrial metabolism and the progression of OA.Additionally, This study sheds light on the involvement of these mitochondrial metabolism-related genes in the pathophysiological processes of OA, offering valuable insights for the future development of targeted therapeutic strategies for this condition.We further analyzed the correlation between these nine mitochondrial metabolism genes and hallmark signaling Enrichment analysis of subtypes.Protein pathway analysis of subtype 1 (A) and subtype 2 (B).GSVA analysis of BP (C) and KEGG (D).
pathways.The results revealed significant associations with various signaling pathways, including unfolded protein response, TNFa signaling via NF-kB, TGF-beta signaling, reactive oxygen species pathway, PI3K-AKT-mTOR signaling, peroxisome, p53 pathway, oxidative phosphorylation, notch signaling, mTORC1 signaling, mitotic spindle, and DNA repair.Mitochondrial unfolded protein response (UPR) is a signaling pathway from mitochondria to the nucleus that is activated to maintain mitochondrial function when misfolded proteins accumulate within the mitochondria.Studies have found that enhanced activation of mitochondrial UPR can reduce chondrocyte death, alleviate joint pain, and lower inflammation levels in synovial fluid (45).TNF-a, a key inflammatory factor in OA pathogenesis, induces chondrocyte death.Research by M J Loṕez-Armada et al. demonstrated that TNF-a and interleukin-1 (IL-1) decrease complex I activity and ATP production, ultimately resulting in mitochondrial depolarization and contributing to cartilage degradation (46).TGF-b3, another growth factor present in OA-affected joints, has shown significant potential in promoting chondrocyte growth and metabolism (47).Specifically, TGF-b3 enhances mitochondrial biogenesis by stimulating mitochondrial fission.TGF-b3-induced mitochondrial fission is mediated by the AMPK signaling pathway.ROS leads to chondrocyte apoptosis, reduced synthesis and activation of matrix metalloproteinases, cartilage degradation, and mitochondrial DNA damage (48).Mitochondria-generated ROS contribute significantly to the pathogenesis of OA.In an OA rat model, artemisinin has been found to activate mitochondrial autophagy by inhibiting the PI3K/AKT/mTOR signaling pathway in IL-1b-induced chondrocytes, leading to the amelioration of disease progression (49).Maintaining a normal number and function of mitochondria relies on the crucial process of mitochondrial biogenesis, which is governed by the central regulator peroxisome proliferator-activated receptor g coactivator-1a (PGC-1a).Sestrin2 (Sesn2), a novel stress-inducible protein, activates AMPK/PGC-1a-mediated mitochondrial biogenesis, thus relieving pain in an OA rat model (50).Curcumin improves OA chondrocytes of OA patients, accompanied by reduced DNA repair capacity (54).The NOTCH signaling pathway is involved in the metabolic processes of cartilage (55).mTORC1 selectively promotes translation of nuclear-encoded mitochondrial mRNA by inhibiting eukaryotic translation initiation factor 4E (eIF4E) binding protein (4E-BP), thereby controlling mitochondrial activity and biogenesis (56).Thus, mitochondria-related genes may play a role in the development of OA through these signaling pathways mentioned above.
Recent research indicates that alterations in inflammation and the immune system contribute to the pathogenesis of OA (57, 58).Moreover, mitochondrial metabolism plays a pivotal role in the metabolism and activation of immune cells.Consequently, we investigated the relationship between mitochondrial metabolismrelated genes and immune cells, identifying significant associations primarily with macrophages, neutrophils, T cells, and dendritic cells.Reprogramming the mitochondrial metabolism of proinflammatory M1 macrophages towards an anti-inflammatory M2 phenotype has emerged as a promising strategy for mitigating OA progression.Through modifications in the mitochondrial metabolism of M1 macrophages, repolarization into M2 macrophages can be achieved, effectively suppressing synovial inflammation and impeding early-stage OA progression (59).Neutrophils are the first immune cells to enter the synovium after joint injury, and their activity is an indispensable prerequisite for the progression of OA (60).Research has revealed that neutrophil elastase (NE) induces chondrocyte apoptosis and promotes the development of OA through the cysteine cathepsin signaling pathway (61).NE, the primary inflammatory protease released by neutrophils, is a key factor in OA.Elevated expression levels of various T cell subtypes, including CD4+ and CD8+ T cells, have been observed in OA patients, indicating their potential role in the pathogenesis of the disease (62).Dendritic cells (DCs) in OA patients also exhibit heightened levels of inflammatory cytokines, suggesting their involvement in the pathogenesis of OA (63).While these immune cell types have been implicated in OA pathogenesis, their precise roles in disease initiation and progression remain poorly understood.Therefore, identifying and comprehending the specific immune cell subtypes involved in OA will provide valuable insights into the disease's etiology and aid in the development of effective treatment strategies.
Inflammation has been shown to play a crucial role in the pathogenesis of OA, particularly through the action of inflammatory cytokines (64).In our study, we identified multiple cytokines that are closely associated with the nine mitochondrial metabolism-related genes.Previous research has reported that interleukin-1b (IL-1b) induces chondrocyte apoptosis by inhibiting SIRT3 expression and mitochondrial autophagy.Mitoquinone (MitoQ5) protects chondrocytes against IL-1binduced oxidative stress and promotes cell survival by upregulating the SIRT3/Parkin-associated autophagy signaling pathway (65).IL-1b increases mitochondrial iron levels to trigger chondrocyte ferroptosis, leading to inflammatory damage (66).In recent years, many studies have demonstrated that alleviating disease progression in OA can be achieved by inhibiting inflammatory cytokines.Artemisinin (AT) has been reported to activate mitochondrial autophagy and alleviate OA by reducing TNFSF11 expression and inhibiting the PI3K/AKT/mTOR signaling pathway (49).Curcumin (CAD) can mitigate inflammation, cartilage degradation, and ferroptosis induced by IL-1b (51).Nodakenin, in an OA mouse model, reduces chondrocyte degradation and inflammatory responses by inhibiting the expression of inflammatory cytokines such as COX-2, IL-1b, and TNF-a (11).Hydrogen sulfide may counteract IL-1b-induced inflammation and mitochondrial dysfunctionrelated cellular apoptosis in chondrocytes by inhibiting the PI3K/ Akt/NF-kB and MAPK signaling pathways (67).
In our study, we acknowledge that there are certain limitations that need to be addressed.Firstly, although we utilized publicly available data for analysis, the reliability of the data could be a potential concern.While efforts were made to ensure data quality and consistency, it is important to note that the original data sources were not generated specifically for our study.Further validation using independent datasets or experimental verification would enhance the robustness of our findings.Secondly, while our pre-modeling approach showed promising results, clinical validation is required to establish the utility and accuracy of the predictive model we developed.Incorporating clinical data and conducting prospective studies would provide valuable insights into the real-world applicability and performance of the model.Lastly, as this study focused on bioinformatics analysis, it is essential to clearly define the directions for future experimental studies.Identifying key genes and pathways associated with mitochondria in osteoarthritis opens up avenues for experimental validation, functional characterization, and mechanistic investigations.Future experimental studies should aim to elucidate the biological relevance of these genes, explore their interactions within cellular processes, and investigate potential therapeutic targets.
In conclusion, the involvement of mitochondrial metabolism in the onset and progression of OA is significant.Targeting mitochondrial dysfunction presents a promising therapeutic approach for OA cases associated with impaired mitochondrial function.However, there is currently a lack of identified drugs that specifically target mitochondrial metabolism in OA.Further research is necessary to comprehend the complex interplay between mitochondrial metabolism and OA.Mitochondrialtargeted therapies hold potential in providing effective preventive and treatment strategies for OA.

1 2 3
FIGURE 1Data processing and enrichment analysis.(A) Data before batch effect removal.(B) Data after batch effect removal.(C) Volcano plot for differential expression analysis.(D) ssGSEA analysis and (E) GSEA analysis reveal a significant decrease in mitochondrial activity in OA. (F) Gene Ontology Biological Process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for mitochondria-related DEGs.(G) WGCNA analysis for mitochondrial scoring.(H) Module Membership (MM) and Gene Significance (GS) values for genes in the black module.(I) Gene Ontology BP and KEGG analysis for genes in the black module.(J) Protein-protein interaction network analysis for genes in the black module.
FIGURE 6 Gene-disease networks and mRNA-miRNA networks.(A) Gene-disease network (red circles represent hub mitochondria-related genes; blue circles represent diseases).(B) Venn diagram displaying the miRNA sets for each hub mitochondria-related DEGs.(C) mRNA-miRNA network (purple shapes represent hub mitochondria-related genes; blue shapes represent miRNAs).

8
FIGURE 8 Identification of molecular clusters and characterization of subtypes based on hub mitochondria-related DEGs in OA. (A) Consensus clustering matrix for k = 2. (B) Cumulative distribution function (CDF) curves for clustering.(C) CDF delta area curves.(D) Tracking plot.(E) Differential analysis of immune cells between subtype 1 and subtype 2 (*p < 0.05, ***p < 0.0001).Drug analysis results for treating patients of subtype 1 (F) and subtype 2 (G).GSEA analysis between subtypes, including BP (H) and KEGG (I) levels.

11
FIGURE 10 WGCNA analysis of subtypes.(A) The optimal soft threshold value is determined to be 4. (B) Merged modules after merging process.(C) Correlation analysis between modules and traits.(D) Heatmap showing the correlation between modules.(E) Scatter plot displaying Module Membership (MM) and Gene Significance (GS) for genes in the grey60 module.(F) BP and KEGG analysis for genes in the module.

12
FIGURE 12Distribution of hub mitochondria-related genes in OA based on single-cell RNA sequencing data.(A) Violin plot showing the distribution of 9 hub mitochondria-related genes in different cell types.(B) The 9 hub mitochondria-related genes produced expression changes in the proposed time series.

TABLE 1
Primers used in this study.